#include <stdint.h>
#include <math.h>
#include "FFT.h"


#define FFT_LENGTH_MAX  8192
#define FFT_PI          3.14159265358979323846f

static float st_sin_table[FFT_LENGTH_MAX];
static float st_cos_table[FFT_LENGTH_MAX];



/*
* Return:      log2(N)
* Parameters:  N:长度
* Description: 求log2(N)
*/
static int32_t iFFTLog2(int32_t N)
{
    int32_t i = 0;

    if(N < 0)
        return 0;

    while(N >>= 1)
        ++i;

    return i;
}

/*
* Return:      void
* Parameters:  ptypeDatas:数据; iLength:长度
* Description: 基2FFT变换
*/
void vFFT(FFTTypeDef *ptypeDatas, int32_t iLength)
{
    static int32_t iLengthLast = 0;
    int32_t m = iFFTLog2(iLength);
    int32_t i, j = iLength / 2, l, iJDown;
    int32_t x, y = iLength / 2, iXNumber, iYNumber;
    FFTTypeDef typeFFTTemp, typeFFTWnk;

    /* 先生成列表，避免反复三角函数计算 */
    if(iLengthLast != iLength)
    {
        for(i = 0; i < iLength; ++i)
        {
            st_sin_table[i] = sin(2.0f * FFT_PI * i / iLength);
            st_cos_table[i] = cos(2.0f * FFT_PI * i / iLength);
        }

        iLengthLast = iLength;
    }


    /* 倒序排序 */
    for(i = 1; i < iLength; ++i)
    {
        if(i < j)
        {
            typeFFTTemp = ptypeDatas[i];
            ptypeDatas[i] = ptypeDatas[j];
            ptypeDatas[j] = typeFFTTemp;
        }
        x = y;
        while(j & x)
        {
            j -= x;
            x >>= 1;
        }
        j |= x;
    }

    /* FFT核心蝶形算法 */
    /* 一重循环对m级的级数做循环 */
    for(l = 0; l < m; ++l)
    {
        /* 当前l级下有iXNumber种类型旋转因子 */
        iXNumber = 1 << l;
        /* 当前l级下有iYNumber组相同的旋转因子 */
        iYNumber = 1 << (m - l - 1);
        /* 每个l级开始做蝶形运算时都是从第0号元素开始 */
        i = 0;

        /* 二重循环计算每个l级下的每个不同旋转因子 */
        for(x = 0; x < iXNumber; ++x)
        {
            /* 参与两个蝶形运算数据的下标 */
            j = i;
            iJDown = j + iXNumber;

            /* 旋转因子变换 */
            typeFFTWnk.real = st_cos_table[x * iYNumber];     //typeFFTWnk.real = cos(2.0 * FFT_PI / N * x * iYNumber);
            typeFFTWnk.imag = -st_sin_table[x * iYNumber];    //typeFFTWnk.imag = -sin(2.0 * FFT_PI / N * x * iYNumber);

            /* 三重循环对每个相应的旋转因子做蝶形运算 */
            for(y = 0; y < iYNumber; ++y)
            {
                /* 蝶形运算 */
                typeFFTTemp.real = ptypeDatas[iJDown].real * typeFFTWnk.real - ptypeDatas[iJDown].imag * typeFFTWnk.imag;
                typeFFTTemp.imag = ptypeDatas[iJDown].real * typeFFTWnk.imag + ptypeDatas[iJDown].imag * typeFFTWnk.real;
                ptypeDatas[iJDown].real = ptypeDatas[j].real - typeFFTTemp.real;
                ptypeDatas[iJDown].imag = ptypeDatas[j].imag - typeFFTTemp.imag;
                ptypeDatas[j].real += typeFFTTemp.real;
                ptypeDatas[j].imag += typeFFTTemp.imag;

                /* 两个蝶形运算的下标 */
                /* 下个相同的旋转因子是在当前号元素标号上加iXNumber的两倍 */
                j += (iXNumber << 1);
                iJDown = j + iXNumber;
            }

            /* 当前l级的下组不同的旋转因子开始标号是在当前旋转因子开始标号上加1 */
            ++i;
        }
    }
}

/*
* Return:      void
* Parameters:  ptypeDatas:数据; iLength:长度
* Description: 基2FFT逆变换
*/
void vIFFT(FFTTypeDef *ptypeDatas, int32_t iLength)
{
    FFTTypeDef *ptypeDataNow = NULL;
    int32_t i;

    /* 变换为共轭对 */
    ptypeDataNow = ptypeDatas;
    for(i = 0; i < iLength; ++i)
    {
        ptypeDataNow->imag = -ptypeDataNow->imag;
        ++ptypeDataNow;
    }

    /* FFT变换 */
    vFFT(ptypeDatas, iLength);

    /* 变换为共轭对后，再除以序列长度 */
    ptypeDataNow = ptypeDatas;
    for(i = 0; i < iLength; ++i)
    {
        ptypeDataNow->real = ptypeDataNow->real / iLength;
        ptypeDataNow->imag = -ptypeDataNow->imag / iLength;
        ++ptypeDataNow;
    }
}
